#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ZFinder.hh"
#include "Rivet/Projections/Beam.hh"

namespace Rivet {
using namespace std;
using namespace Cuts;

class PP_2 : public Analysis {
private:
	vector<Histo1DPtr> _h_dists;
	vector<Histo1DPtr> _h_histos;
	Scatter2DPtr _h_test;

public:
	PP_2() : Analysis("PP_2")
	{
	}

	void init()
	{
		Beam beams;
		FinalState fs;
		ZFinder zfinder(fs, Cuts::open(), PID::ELECTRON, 60.0 * GeV,
		    100.0 * GeV, 0.0 * GeV, ZFinder::CLUSTERNODECAY);

		addProjection(zfinder, "ZFinder");
		addProjection(beams, "Beams");

		_h_dists.resize(13, Histo1DPtr());
		_h_dists[0] = bookHisto1D("A0", 50, 0, 600);
		_h_dists[1] = bookHisto1D("A1", 50, 0, 600);
		_h_dists[2] = bookHisto1D("A2", 50, 0, 600);
		_h_dists[3] = bookHisto1D("A3", 50, 0, 600);
		_h_dists[4] = bookHisto1D("A4", 50, 0, 600);
		_h_dists[5] = bookHisto1D("A5", 50, 0, 600);
		_h_dists[6] = bookHisto1D("A6", 50, 0, 600);
		_h_dists[7] = bookHisto1D("A7", 50, 0, 600);
		_h_dists[8] = bookHisto1D("fL", 50, 0, 600);
		_h_dists[9] = bookHisto1D("fR", 50, 0, 600);
		_h_dists[10] = bookHisto1D("f0", 50, 0, 600);
		_h_dists[11] = bookHisto1D("ptz", 50, 0, 600);
		_h_dists[12] = bookHisto1D("divide", 50, 0, 600);

		/*
		_h_dists[1] = bookHisto1D("A1", logspace(100, 1.0,
		    0.05 * sqrtS()));
		_h_dists[2] = bookHisto1D("A2", logspace(100, 1.0,
		    0.05 * sqrtS()));
		_h_dists[3] = bookHisto1D("A3", logspace(100, 1.0,
		    0.05 * sqrtS()));
		_h_dists[4] = bookHisto1D("A4", logspace(100, 1.0,
		    0.05 * sqrtS()));
		_h_dists[5] = bookHisto1D("A5", logspace(100, 1.0,
		    0.05 * sqrtS()));
		_h_dists[6] = bookHisto1D("A6", logspace(100, 1.0,
		    0.05 * sqrtS()));
		_h_dists[7] = bookHisto1D("A7", logspace(100, 1.0,
		    0.05 * sqrtS()));
		_h_dists[8] = bookHisto1D("fL", logspace(100, 1.0,
		    0.05 * sqrtS()));
		_h_dists[9] = bookHisto1D("fR", logspace(100, 1.0,
		    0.05 * sqrtS()));
		_h_dists[10] = bookHisto1D("f0", logspace(100, 1.0,
		    0.05 * sqrtS()));
		_h_dists[11] = bookHisto1D("ptz", logspace(100, 1.0,
		    0.05 * sqrtS()));
		_h_dists[12] = bookHisto1D("divide", logspace(100, 1.0,
		    0.05 * sqrtS()));
		*/

		_h_histos.resize(4 , Histo1DPtr());
		_h_histos[0] = bookHisto1D("thetastar", 100, -1.0, 1.0);
		_h_histos[1] = bookHisto1D("phistar", 90, 0.0, 360.0);
		_h_histos[2] = bookHisto1D("thetastar_ptw20", 100, -1.0, 1.0);
		_h_histos[3] = bookHisto1D("phistar_ptw20", 90, 0.0, 360.0);

		_h_test = bookScatter2D("test");
	}

	void analyze(const Event& event) {
		const double weight = event.weight();
		const ZFinder& zfinder = applyProjection<ZFinder>(event,
		    "ZFinder");

		if (zfinder.bosons().size() != 1)
			vetoEvent;

		const ParticlePair& beams = applyProjection<Beam>(event,
		    "Beams").beams();
		FourMomentum pb1(beams.second.momentum()),
		    pb2(beams.first.momentum());
		FourMomentum lm = zfinder.constituents()[1].momentum();

		FourMomentum pl(lm);
		FourMomentum pll(zfinder.bosons()[0].momentum());
		LorentzTransform cms(-pll.boostVector());
		Matrix3 zrot(pll.p3(), Vector3(0.0, 0.0, 1.0));

		pl = cms.transform(pl);
		pb1 = cms.transform(pb1);
		pb2 = cms.transform(pb2);

		Vector3 pl3 = pl.p3();
		Vector3 pb13 = pb1.p3();
		Vector3 pb23 = pb2.p3();
		pl3 = zrot*pl3;
		pb13 = zrot*pb13;
		pb23 = zrot*pb23;
		Vector3 xref(cos(pb13.theta()) > cos(pb23.theta()) ?
		    pb13 : pb23);
		Matrix3 xrot(Vector3(xref.x(), xref.y(), 0.0),
		    Vector3(1.0, 0.0, 0.0));
		pl3 = xrot * pl3;

		double ptz(zfinder.bosons()[0].pT() / GeV);
		double thetas(pl3.theta()), phis(pl3.phi());
		double costhetas(cos(thetas)), sinthetas(sin(thetas));
		double cosphis(cos(phis)), sinphis(sin(phis));

		if (phis < 0.0)
			phis += 2.0 * M_PI;

		_h_histos[0]->fill(costhetas, weight);
		_h_histos[1]->fill(phis * 180.0 / M_PI, weight);

		if (ptz > 20.0) {
			_h_histos[2]->fill(costhetas, weight);
			_h_histos[3]->fill(phis * 180.0 / M_PI, weight);
		}

		_h_dists[0]->fill(ptz, (10.0 / 3.0 * (1.0 - 3.0 *
		    sqr(costhetas)) + 2.0 / 3.0)* weight);
		_h_dists[1]->fill(ptz, 10.0 * sinthetas * costhetas * cosphis *
		    weight);
		_h_dists[2]->fill(ptz, (10.0 * sqr(sinthetas) * (sqr(cosphis) -
		    sqr(sinphis))) * weight);
		_h_dists[3]->fill(ptz, 4.0 * sinthetas * cosphis* weight);
		_h_dists[4]->fill(ptz, 4.0 * costhetas* weight);
		_h_dists[5]->fill(ptz, 4.0 * sinthetas * sinphis* weight);
		_h_dists[6]->fill(ptz, 10.0 * costhetas * sinthetas * sinphis*
		    weight);
		_h_dists[7]->fill(ptz, 10.0 * sqr(sinthetas) * cosphis *
		    sinphis* weight);
		_h_dists[8]->fill(ptz, (0.5 * sqr(1.0 - costhetas) - (1.0 - 2.0 *
		    sqr(costhetas)))* weight);
		_h_dists[9]->fill(ptz, (0.5 * sqr(1.0 + costhetas) - (1.0 - 2.0 *
		    sqr(costhetas)))* weight);
		_h_dists[10]->fill(ptz, (5.0 * sqr(sinthetas) - 3.0)* weight);
		_h_dists[11]->fill(ptz, 1* weight);
	}


	void finalize() {
		foreach (Histo1DPtr histo, _h_histos) {
			scale(histo, crossSectionPerEvent());
		}

		divide(_h_dists[2], _h_dists[11], _h_test);

		/*
		   std::vector<YODA::ProfileBin1D>::iterator it0, it1;
		   std::vector<YODA::ProfileBin1D> bins0, bins1;
		   int i;
		   bins0 = _h_dists[0][2]->bins();
		   bins1 = _h_dists[0][11]->bins();
		   double y;
		   for (it0 = bins0.begin(), it1 = bins1.begin(), i = 0; it0 != bins0.end() && it1 != bins1.end(); i++, it0++, it1++) {
		   if ((*it1).sumWY() == 0)
		   y = 1;
		   else
		   y = (*it0).sumWY() / (*it1).sumWY();
		   _h_dists[0][12]->fillBin(i, y, 1.0);
		   }
		   */
	}
};

DECLARE_RIVET_PLUGIN(PP_2);
}
